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ABSTRACT 

A nonlinear kinetic theory, combining cosmic-ray (CR) acceleration in supernova remnants (SNRs) 
with their gas dynamics, is used to re-examine the nonthermal properties of the remnant of SN 1987A 
for an extended evolutionary period of 5-50 yr. This spherically symmetric model is approximately 
applied to the different features of the SNR which consist of (i) a blue supergiant wind and bubble, 
and (ii) of the swept-up red supergiant (RSG) wind structures in the form of an H II region, an 
equatorial ring (ER), and an hourglass region. The RSG wind involves a mass loss rate that decreases 
significantly with elevation above and below the equatorial plane. The model adapts recent three- 
dimensional hydrodynamical simulations by Potter et al. in 2014 that use a significantly smaller 
ionized mass of the ER than assumed in the earlier studies by the present authors. The SNR shock 
recently swept up the ER, which is the densest region in the immediate circumstellar environment. 
Therefore, the expected gamma-ray energy flux density at TeV energies in the current epoch has 
already reached its maximal value of 10“^^ erg cm“^ s“^. This flux should decrease by a factor of 
about two over the next 10 years. 

Subject headings: acceleration of particles — ISM: individual objects (SN 1987A) — ISM: supernova 
remnants — X-rays: individual (SN 1987A) — gamma rays: ISM 


1. INTRODUCTION 

The supernova (SN) that occurred in 1987 in the 
nearby Large Magellanic Gloud was the hrst object of its 
kind whose evolution in the radio to X-ray range has been 
resolved as a function of time. The study of SN 1987A 
continues and includes, in particular, extensive observa¬ 
tions in__verYAiigh__energy_XVH^E > 100 GeV) 7 -rays 
(e.g. IH.E.S.S. Collaborationi[2M5ll . 

The present work is a critical re-examination and ex¬ 
tensi on of the studies bviBer ezh ko fc KsenofontovI (I2000L 
l200(ifl and iBerezhko et al.llj201ll hreferred to as BKVll 
in the following) concerning the properties of the non¬ 
thermal emission of SN 1987A. As in those previous in¬ 
vestigations, our framework is a nonlinear kinetic theory 
of cosmic-ray (CR) acceleration in SN remnants (SNRs). 
This theory couples the particle accelerat ion process with 
the h ydrodynamics of the thermal gas (jBerezhko et al.l 
199(3) and connects it w ith the gamma-ray emission (e.g. 
Berezhko fc VblljI2000D. The application of this theory 

to ind ividual SNRs ('see lVolkl2004lBerezhkol[2005Ll200^ 

12014 for reviews) has demonstrated its ability to ex¬ 
plain the observed SNR properties and radiation spec¬ 
tra. Combining the theoretical model with the observed 
synchrotron spectra predicts new effects like the large 
degree of magnetic field amplification that leads to the 
observed concentration of the highest-energy electrons in 
a very thin shell just behind the forward shock into the 
circumstellar medium (CSM). 

The evolution of radio and X-ray emission at earlier 
times, also implied in the present paper, has been de¬ 
scribed by, e.g. BKVll, who showed, in particular, that 
the evidence for strong shock modification comes primar- 
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ily from radio data. 

Efficient acceleration of the CR proton component is 
needed to produce significant shock modification lead¬ 
ing to a soft and concave CR electron spectrum in SN 
1987A, which well fits the observed nonthermal radio 
emission and X-ray spectra if the downstream magnetic 
field strength is as high as ~ 15 mG. Such a high 
field leads to significant electron synchrotron losses that 
cut off the high-frequency X-ray part of the synchrotron 
spectrum, consistent with X-ray observations. The re¬ 
quired magnetic held strength can presumably be at¬ 
tributed to its nonlinear amplihcation ne ar the SN R 
shock by the GR acceleration process itself (lBellll2004) . 

As a consequence of the efficient production of the 
CR nuclear component, which is accompanied by strong 
magnetic held amplihcation, SN 1987A is expected to 
be a potential source of y-rays for the H.E.S.S. instru¬ 
ment. However, the hux of TeV emission predicted in 
BKVll exceeds the upper li mit obtained by H.E.S.S. 
(jH.E.S.S. Collaborationlf20T^ . It is suggested here that 
this inconsistency is the result of an overestimate of the 
CSM gas density, in particular, of the mass of the equa¬ 
torial ring (ER). 

The present work uses the same canonical values (e.g. 
lMcCravlll993l) for the stellar ejecta mass, Mej, distance, 
d, hydrodynamic explosion energy, Egn, and ejecta ve¬ 
locity distribution as in BKVll, but takes int o account 
the de tailed radio continuum obs ervations by iNg et al.1 
(|2013l) and iZanardo et al.l (|20l4) . In particular, the 
present work makes use of recent extensive CSM mod¬ 
eling in ter ms of three-dimensio nal hydrodynamical sim¬ 
ulations by iPotter et al.1 (|20l4 ) ( in the sequel referred 
to as Potter 14). An approximate prediction of the fu¬ 
ture nonthermal emission from SN 1987A for the coming 
decades is given. The densest part of the CSM, which lies 
within the range |0| < 20° of the elevation angle 0 around 
the equatorial plane, is then considerably less dense than 
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what was adopted in BKVll. This is the main reason for 
the re-examination of the CR production in SN 1987A 
and the associated nonthermal emission. 

2. SN 1987A AND ITS CIRCUMSTELLAR ENVIRONMENT 

To study the propagation of the SN shock through 
the CSM, the results of Potterld for the angular range 
|0| < 20° relative to the equatorial plane are usetfl. 
The most efficient CR and nonthermal emission produc¬ 
tion presumably takes place within this regi on. This 
is rou ghly consistent with radio observations (iNg et al.l 
1201311 . The adopted radial profile of the gas number 
density Ng = p/nip in this region is represented in 
Figure [TJ Within the selected elevation range, it con¬ 
sists of several different morphological structures: (i) the 
wind bubble of the blue supergiant (BSG) progenitor star 
(|Chevalier fc FranssonlllQS^ at r < Rq = 4.5 x 10^^ cm 
with gas density Vp- = 0.29 c m~^, (ii) the H II region 
(jChcvalier fc Dwarkadasi Il995f l at i?c < r < Rhg = 
8 X 10^^ cm with Ng = 280 cm“^, (iii) the so-called 
hourglass region at Rhg < r < i?w = 1-5 x 10^® cm with 
Ng = 10 cm“®, and (iv) the free red supergiant (RSG) 
wind region at R > i?w with Ng = 10(r/Rw)^ cm“®; the 
properties of structures (iii) and (iv) directly follow from 
Potterl4. Within the smaller elevation angle region of 
|0| < 4.5°, the same radial profile includes the equatorial 
ring inside the H II region (see Figure [T]). Its gas number 
density is chosen here to be distributed according to the 
relation 

iVg = iVgmexp[-(r-RER)V^ER]: (1) 

where Ngm ~ MER/(47r®/^mpi?|E^ER) is the central 
(maximal) density of the ER, and Mer, ReRj and Zer 
denote the total mass, radius, and width of the ER, re¬ 
spectively. Below, the values Mer = O.O58M0, Rer = 
6.4 X 10^^ cm, and /er = 0.12Rer are used. These pa¬ 
rameter values are taken in order to fit the observed 
shock and radio-emission dynamics. They turn out to 
be consistent with the values used by Potterl4. 

The CSM described above differs significantly from the 
CSM adopted in BKVll. First of all, the ER mass, al- 
beit that of the ioni zed material only, Mer = O.O58M0 
(|Mattila et al.l [201^ . is considerably smaller than that 
used by BKVll (Mer = O. 5 M 0 ). Second, the CSM 
behind the ER is less dense. These factors, as demon¬ 
strated below, lead to a considerable reduction of the ex¬ 
pected nonthermal emission, in particular, of the 7 -ray 
emission, compared with the results of BKVll. 

3. PARTICLE ACCELERATION MODEL 
3.1. Shock Approximation 

The propagation of the forward SN shock through the 
CSM is modeled in the spirit of a spherically symmetri¬ 
cal approach. It approximates the shock and its effects 
as the weighted sum of two independent spherically sym¬ 
metric shocks, propagating into, respectively, two differ¬ 
ent radial gas number density profiles of the CSM, shown 
in Figure [TJ the first profile (region 1 ) belongs to an az- 
imuthally symmetric region of elevation angles 9 near the 

^ The equatorial-to- polar density ratio is 20:1 
(IBlondin fc LundQvistlll993l) . 



Fig. 1.— Radial profiles of CSM gas density Ng. The partly 
overlapping thick dashed and solid lines, respectively, correspond 
to region 1 (4.5° < |6| < 20°) and region 2 (|6| < 4.5°) in elevation 
angle 6 . The radial extent of the various morphological structures, 
summarized in section 2, is indicated by the vertical dotted lines. 

equatorial plane 4.5° < |0| < 20°, and the second, anal¬ 
ogous profile (region 2 ), corresponds to the innermost 
|0| < 4.5°. Then, each quantity Q, characterizing the 
number of accelerated CRs and the amount of emission 
produced by these CRs, is determined by the relation 

Q = Qi(/i —/2) + Q2/2, ( 2 ) 

where Qi ,2 are the spherically symmetric values corre¬ 
sponding to the profiles that characterize the two regions 
1 and 2, respectively, and /i = 0.34 and /2 = 0.08 denote 
the filling factors in the solid angle of these regions. 

CR production by the SN shock at high latitudes 
| 0 | > 20 ° is neglected here because the gas density in 
this region is considerably lower (see Potterl4). We also 
neglect CR production at the reverse shock for the rea¬ 
sons given in BKVll. 

The question is, of course, to what extent such an ap¬ 
proximate treatment remains consistent as a function of 
time with the real SNR shock sweeping across the real 
CSM to reach larger and larger radial distances. A pri¬ 
ori the decrease of the gas density with elevation angle, 
symmetric to the equatorial plane, suggests this to be 
a roughly stable process. However, effects like the en- 
gulfment of the ER clearly imply some non-radial shock 
propagation aspects. Such effects are also apparent in the 
three-dimensional simulations of Potterl4. In addition , 
the ring is clumpy, even though iFransson et al.l (1201511 
found indications that these hot spots are now gradually 
dissolving. This type of effect should primarily influence 
the detailed time dependence of the particle acceleration 
and hadronic 7 -ray emission rather than the global accel¬ 
eration properties of the systenfl. Therefore, in a gross 
sense, the used mosaic of spherical shocks appeears to be 
an adequate overall approximation. 

3.2. Field Amplification, Injection Rate, and Electron 
to Proton Ratio 

As in BKVll, the magnetic field strength, Bq, given 
by the expression 

Bo = v'27r X lO-VoK^ (3) 

F or di scussions of such effects, h owever, see IBerezhko et al.l 
1120131) and IGabici &: AharonianI 1120141) . 
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Fig. 2.— Shock radius Rb and shock speed 14 (solid lines), gas 
number density Ng (dash-dotted line) and upstream magnetic field 
Bo (dashed line) at the current shock position as a function of time 
since SN explosion, for region 2. The dotted vertical line marks 
the current epoch. The observed radius Rs an d speed 14 of th e SN 
shock, as determined by radio observations lINg et al.l 1201311 are 
shown as well. The scaling values are R[ = R^ = 3.1 X cm 

and Vi = 28000 km s“^. 

is used. Here, denotes the SN shock speed and Bq is 
the field far upstream, presumably amplified by the CRs 
of the highest energy. In the same sense, po is the mass 
density far upstream. The high downstream magnetic 
field Bd = Bq X a Ki 10 mG, where a denotes the total 
shock compression ratio, is req uired to reproduce the ob¬ 
serve d radio and X-ray spectra fBerezhko fc KsenofontovI 
I2006D . The calculation of the s hock radius and speed 
K ag ain follows the scheme of (iBerezhko fc KsenofontovI 
I2006D , but see also BKV11, as does the evaluation of the 
proton injection rate T^iVg and of the electron-to-proton 
ratio ibep- 


4. RESULTS AND DISCUSSION 

Eigure[2]shows i?s and 14, as the shock propagates in 
the GSM corresp onding to regio n 2, together with the 
latest radio data (iNg et al.ll2013D . In the case of region 
I, the shock speed time profile 14(t) does not contain the 
local minimum around t = 8000 days; this is the main 
difference from the results presented in Figure [2j Itera¬ 
tively fitting the theoretical quantities r]{t) and Kep{t) to 
the spatially i ntegrated radio synchrotron spectra up to 
the year 2013 (INg et al.ll2013D leads to a constant value 
for Kep{t) = 3 X 10“^, whereas the value r]{t) r:! 3 x 10“^ 
at t Ri 26 year (Figure [3^) is due to the assumption that, 
leaving the H B region, the nuclear injection fraction 
should, after t Ri 30 year, go back to its value before the 
age of 10 years. During the ages between about 10 and 
30 years the compressed and largely azimuthal magnetic 
field of the H B region should have depressed nuclear 
injection, consistent with the softening of the spatially 
integrated radio spectrum (with index a{t)) as shown in 
Figure [3 )d. This figure also shows that outside the H B 
region, the shock is significantly modified relative to a 
pure gas shock for which the total shock compression ra¬ 
tio tr and the subshock compression ratio CTs would both 
have a value of 4. 

The adopted value of ri{t) determines the amount of 
shock modification, in particular, the decrease of the sub¬ 
shock compression ratio crs(t) relative its value of 4 for 
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Fig. 3. — (a) Shock compression ratio a and subshock compres¬ 
sion ratio (Tg (solid lines), proton injection fraction rj (dashed line), 
and electron-to-proton ratio Kep (dash-dotted line) as functions of 
time, (b) Self-consistent spectral index a of the integrated nonther- 
mal radio emissio n as a function of time together with observational 
data from ATCA llZanardo et a n[2oTq) and the latest combination 
of ATCA and ALMA llZanardo et al.ll2014l) . 

an unmodified shock. Since CRs with relatively small 
energies are produced near the subshock, the value CTs(t) 
directly determines the shape of the electron energy spec¬ 
trum at energies below ri 1 GeV — which produce syn¬ 
chrotron emission in the radio range — and vice versa. 
Therefore, as in all similar cases, the proton injection 
rate r]{t) is inferred from a fit of the resulting theoretical 
spectral index a(t) of the integrated radio synchrotron 
emission to the observed value. The quality of this fit 
can be ascertained from Figure [Sb, where the best cal¬ 
culated time profile a{t) is presented together with the 
observational data. 

From the epoch t > 6 year onwards th e value of a(t) 
decre ases until the epoch t ss 23 year (jZanardo et al.l 
l2?Bnl) and is assumed to start to increase again at t > 23 
year. This a ppears to be consi s tent w ith recent mea- 
surements bv lindebetouw et al.l (|20I4D and, especially, 
iZanardo et al.l ( 20141) . However, the increase is not as 
fast as one would expect from the behavior of the local 
This occurs because at subsequent epochs, the gas 
number density and, consequently, also the amount of 
freshly injected electrons and the magnetic field strength, 
are much smaller than when the shock was in the H B re¬ 
gion. Therefore, the synchrotron emission of those latter 
electrons, convected downstream in that higher magnetic 
field, remains dominant. For similar reasons, the varia¬ 
tion in the shock compression ratio a at an age 10,000- 
12,000 days, caused by the rapid decrease of the upstream 
gas number density Ng, does not significantly affect the 
integrated radio synchrotron emission flux, and the spec¬ 
tral index a remains essentially independent of a. 

Using these educated guesses for ry and ibep tbe calcu¬ 
lated flux of radio emission Sn at frequency v = 0 GHz 
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Fig. 4. — Calculated flux of radio emission at frequency v = 
9 GHz as a fu nction of time, together with the ATCA data 
(INg et aLII2013l) . Dashed and solid lines represent the contribu¬ 
tion of region 1 and of the sum from regions 1 and 2 , respectively. 

is presented in Figured together with the obs ervational 
data obtained with the ATCA instrument (|Ng et alJ 
1201311 . For those times where the radio flux has also 
been measured, the calculation is perfectly consistent 
with the observations. According to the calculation, the 
rapid growth of radio emission in the epochs t < 30 year 
is due to the increase of the number of accelerated CR 
electrons, which is proportional to the swept-up mass 
within the H II region. For the adopted value of the 
outer boundary of the H II region, Rhg = 8 x 10^^ cm, 
the SNR shock reaches this boundary after t ss 30 year, 
and the peak of radio emission is reached. Even at later 
epochs, t = 30 —40 years, the radio synchrotron emission 
will still be dominated by the contribution of the swept- 
up matter of the H II region. The flux of radio emission 
is expected to decrease gradually by a factor of about 10 
during the 10 years after 2017. Due to its considerably 
lower gas density, the contribution of the swept-up hour¬ 
glass matter will become essential only for t > 43 year, 
when the radio emission starts to grow slowly again. 

The calculated y-ray energy flux density above 3 TeV 
as a function of time, shown in Figure El is dominated 
by the 7 r*^-decay component at all energies. Since a sig- 
nihcant part of the shock surface is expected to be tan¬ 
gential, and therefore to not efficiently inject/accelerate 
nuclear CRs, the overall number of accelerating CRs is 
no rmalized by a facto r of /re = 0 . 2 , as previoously argued 
bv iVbIk et all ([2n?ill and BKVll. 

As is clear from Figure El the region 1 , which con¬ 
tains the ER, the densest structure, contributes domi¬ 
nantly only during the shock propagation through the 
ER, which is during days 7000-10,000. 

According to Figure El the maximal energy flux den¬ 
sity of TeV emission Ri 7 x 10“^ eV cm“^ s“^ was 
achieved at day 9000, and after that epoch it decreased 
continuously due to the decrease of the CSM gas density. 
Since the radial gas density profile is much thinner com¬ 
pared with the one used earlier BKVll, the peak value 
of the expected flux of TeV emission is now lower and 
the TeV emission expected for the future is considerably 
lower. 

The most recent upper limit for the TeV emission 
obtained by the H.E.S.S. telescopes during the period 
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Fig. 5. — Integral 7 -ray energy flux density above 3 TeV 
from SN 1987A as a function of time. The H.E.S.S. 
(|H.E.S.S. Collaboration! I2Q15I) upper limit, corresponding to the 
observational period 2005-2012 year, is shown as well. 

2005-2012 (IH.E.S.S. CollaborationI (see Figure El 
is roughly consistent with this prediction. 

According to the calculation (see Figure El, the most 
promising time for the detection of SN 1987A in TeV 7 - 
rays is the 10 year period from 2008 to 2018. At later 
epochs, SN 1987A should be detectable in the VHE range 
only by an instrument with a higher sensitivity than that 
of H.E.S.S. 

This work has been supported in part by the Rus¬ 
sian Foundation for Basic Research (grants 13-02-00943 
and 13-02-12036) and by the Council of the President 
of the Russian Federation for Support of Young Scien¬ 
tists and Leading Scientific Schools (project No. NSh- 
3269.2014.2). 

5. APPENDIX 



Fig. 6. — Spatially integrated synchrotron spectral energy flux 
density of SN 1987 A, calculated for five evolu tionary epochs. The 
ATCA radio data JZanardo et al.ll2010l . 1^1411 for four epochs are 
shown as well, together with the Chandra X-ray flux JPark et al.l 
1200711 data for two epochs (both in the soft energy range at = 10^^ 
Hz and, connected by dotted lines, in the hard energy range at 
u = 6x 10^'^ Hz). 
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Fig. 7.— Spatially integrated 7 -ray spectral energy flux den- 
sity, calculated for five epoch s. The recent H.E.S.S. upper limit 
dH.E.S.S. CollabQratiQnll2015l) is shown as well. 
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